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By carrying out density functional theory (DFT) analysis within the Keldysh non-equilibrium 
Green's function (NEGF) framework, we investigate the quantum transport properties of Au-Ceo- 
' An molecular junctions from the first principles. We briefly review the NEGF-DFT formalism and 

\ present some of our data. We found that at small electrode separations (~ 10 — 13 A), the Au- 

■ Ceo-Au junctions show metallic behaviour with G ~ (2e^//i) and / ~ 1 — 3/iA in the linear regime 
^ ' (|Vbias| < 0.3V). The physical mechanism is the resonance tunnelling through partially occupied 
^ , states originated from the lowest unoccupied molecular orbital (LUMO) of the free Ceo. We also 

■ found that the charge transfer from the Au electrodes to the Ceo molecule can be controlled by gate 
_ ' potential. 

\o ■ 

(N : 

I I. INTRODUCTION 

^ : 

^ . Using molecules as functional units for electronic device application [1] is an interesting perspective and a possible 
' goal of nano-electronics. Work in this field has clearly demonstrated that many of the important molecular device 
-^-^ , characteristics relate specifically to a strong coupling between the atomic and the electronic degrees of freedom, for 
^ • a popular introduction see reference [2] . However, from a theoretical point of view, the accurate prediction of the 
I t ' properties of atomic and molecular scale devices - including the true I-V curves with as few adjustable parameters as 
d • possible - still represents a formidable challenge [3-5] despite the advances and wide-spread application of large scale 
^ ' ah initio modeling based on the density functional theory (DFT) over the last two decades 

I . Recall that most of the previous DFT-based ab initio simulations [6,7] solve only two kinds of problems: (i) finite 
systems such as isolated molecules and clusters, as in quantum chemistry; (ii) periodic systems consisting of super- 
cells, as in solid state physics. However, a molecular electronic device is neither finite nor periodic. Typically, it has 
^ open boundaries which are connected to long and different electrodes extending to electron reservoirs far away, and 
the external bias potentials are applied to these reservoirs. In other words, calculations of finite or periodic systems 
do not include the correct boundary conditions for a nano-scale device in the quantum transport regime. Therefore, 
a new formalism for electronic analysis is required to carry out the first principle transport modeling for molecular 
electronics. 

In this article, we briefly outline a formalism that combines the density functional theory with the Keldysh non- 
\^ , equilibrium Green's functions (NEGF) so that non-equilibrium properties of the quantum transport regime can be 
0^ ' predicted from atomistic approach without any phenomenological parameters [8,9]. To show how it works, we apply 
, our NEGF-DFT approach to investigate a Ceo molecule connected to two Au electrodes. The rest of the article is 

■ organized as follows: in section II, we present the NEGF-DFT formalism and, in section III, we present the transport 
properties of the Au-Ceo-Au tunnel junctions and compare them with the properties of the Al-Cgo-Al ones. Section 

. IV is reserved for a short summary. 
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FIG. 1. Schematic plot of an Au-Ceo-Au molecular tunnel junction. The Au electrodes consist of repeated unit cells extending 
to 2 = ±00 of the horizontal axis z and its surface is oriented along the (100) direction. The lower panel shows the calculated 
charge density at equilibrium (/i; — /ir)- Note that perfect matching is obtained across the boundaries between the leads and 
the central scattering region. 



II. NEGF-DFT FORMALISM 



The simplest model of a molecular device is schematically shown in figure 1 where a Cgo molecule is contacted by 
two semi-infinite Au electrodes which extend to z = ±00 where bias voltage is applied. To calculate the electronic 
states of such a device, two problems should be solved. First, the infinitely large problem (due to the electrodes) must 
be reduced to something manageable on a computer, i.e., one has to solve an open boundary problem. Second, one 
has to find the charge density p{r) of the molecule and electrodes provided there is a bias voltage across the device, 
and this is a non-equilibrium problem. From the computational point of view, it is convenient to divide our system 
into three parts: a scattering region with some portion of electrodes and two metal electrodes extending to z = ±00, 
as shown by the vertical lines in figure 1. 

To reduce the infinitely large problem to that defined inside the scattering region (see figure 1), we notice that the 
effective Kohn-Sham (KS) potential Vi2ff[p(r)] deep inside the left or right lead is very close to the corresponding bulk 
KS potential. This fact makes the boundary conditions to be written in the following form [8] : 

{Vi,off'(r) = V^,buik(r), z < zi, 
Vc,off'(r), zi < z < Zr, (1) 

off (r) = V;., bulk (r), z>Zr, 

where the planes z = z^^) are the left (right) limits of the scattering region (see figure 1), and V;^buik(r) and buik(r) 
are known. In practice, within the DFT, we only need to match the Hartree potential Vfj[p(r)] at the boundaries. 
This can be accomplished by solving the Poisson equation on C/jj(r) by use of a 3D real space grid and boundary 
conditions at z — z^^), and Vfj(r) = e?7jj(r). Once the Poisson equation is solved, Vfj(r) is matched, and by use of 
the DFT we can show that Vcs{r) is also matched. Thus, equation (1) is valid. 

In addition, the 3D real space Poisson solver allows us to deal with any gate potentials, Vg. In our model, they 
provide just additional boundary conditions for the electrostatics, so they change C/jj(r). Furthermore, once 14fr(r) is 
matched across the boundaries, the charge density p(r) is automatically matched at the boundaries. The lower panel 
of figure 1 shows the charge density: although the densities in the leads and in the scattering region are calculated 
separately, they match perfectly across the boundaries. 
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The above "screening approximation" [8] allows us to reduce the infinitely large Hamiltonian of the system to 
the effective finite one defined inside the scattering region. In this approximation, we neglected any infiuence the 
scattering region may give to the leads, but if the portion of the leads included inside the scattering region is long 
enough, such an approximation is well controlled. On the other hand, the semi-infinite leads do contribute to the 
scattering region, and this contribution is included through the self-energies S;(r)(ii^) in the Green's function G{E) 
of the scattering region, for example, 

CR{E) = {E-Ho- Fps -V^- Kc - Sf - S«) , (2) 



V^ = eU^, A[/H(r) = -47rp(r). (3) 

Here, Hq is the kinetic energy of a valent electron, and the atomic cores of the scattering region are fixed in space and 
described by pseudo-potentials, so we have V^ps, and ]4ff[p] = Vfj[p] -I- V^c[p] is the effective KS potential obtained 
within the DFT. For example, the exchange correlation potential V^c is taken in a simple local form following the 
polynomial formula suggested in reference [10], and the (retarded) self-energies due to the semi-infinite left and right 
electrodes are calculated following an iterative technique described in references [8,11]. 

However, to calculate the charge density p{r) away from equilibrium due to the bias voltage 14, we have to apply 
the NEGF method [12]. The density matrix of the scattering region is calculated from the so-called lesser Green's 
function G^{E) of the scattering region as follows: 

p=-^ [dEG<{E), where G<{E) =G^{E)i:<{E)G'^{E). (4) 
27^^ J 

Here, G^^"^\E) is the retarded (advanced) Green's function of the scattering region. Note that G^{E) is defined 
through the Keldysh equation [12]. The lesser self-energy S<[/;, f^] can be evaluated within the mean field theory as 
follows: 



Y.<{E)^ifi{E, fii)riiE)+ifriE, fir)rr{E), (5) 

l{r) = n^Kr) - ^l\r)) Sf^,*,^ 



where F; (^-j = i{Yif,,^ — Sw^-,) and Sw'.^^' is the retarded (advanced) self-energy of the left (right) electrode, and 



fi(^r){E) are the corresponding Fermi distribution functions. Note that fr] is more than a simple Fermi 

distribution: a fact reflecting the non-equilibrium nature of the problem, jii ^ jir- 

In our numerical procedure, we use a LCAO minimal s, p, d atomic basis set (a real space fireball linear combination 
of atomic orbitals [7]) to expand the electron wavefunction and define the effective Hamiltonian matrix F: 

Ho + Vps + V^ + ^ Fjk [ p{r) ] 

Here, to calculate Vpsj^, the atomic cores are described by standard norm conserving nonlocal pseudo-potentials [13]. 
We evaluate G^ by direct matrix inversion, G"^{E) = G^^E), and we construct p{r) from numerical evaluation of 
the integral in equation (4) using a suitable contour in the complex plane E [8]. Once p{r) is obtained, we evaluate 
the KS potential and iterate the above procedure, equations (2), (3), and (4), untill numerical convergence is reached. 

Thus, we construct the self-consistent p{r) by the NEGF technique (it takes care of the non-equilibrium statistics) 
and we calculate the Hartree potential by solving the 3D Poission equation directly (equation (3) includes all the 
external fields as the electrostatic boundary conditions). As a result, our NEGF-DFT formalism is able to solve 
charge transport problems in addition to the conventional electronic structure calculations: we know both the density 
of states DOS{E) and the density matrix G<(r, r', E)/2m of the scattering region. Moreover, we can easily calculate 
the electron charge on the molecule at nozero Vf, and Vg. Of course, the density matrix p is normalized to the total 
number of valent electrons in the scattering region, 

Tr(/56') = TVe = const, 

where S is the overlap matrix of the atomic orbitals we use. Then, one can introduce the (average) number of electrons 
in the molecule by calculating this trace over the molecular orbitals belonging to the molecule only: 



Trmoi(p5) = Are,moi(V6, Vg). 
;,moi(H, Vg) and NfXoX of 

charge on the molecule in the units of e. 



The differcncx! between the calculated A^e^moKH, Vg) and Nf']^^-^ of a free molecule can called an (average) excess 



g = A^e,mol(H, V,)-7Vi°L- 
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III. AU-Cso-AU MOLECULAR TUNNEL JUNCTION 

In this section we report our analysis on the transport properties of an Au-Ceo-Au molecular tunnel junction 

calculated by use of our NEGF-DFT electronic package as discussed above. The device structure is shown in the 
upper panel of figure 1. So far a considerable amount of theoretical and experimental effort has been devoted to 
investigate transport properties of Ceo and other fullerene molecules [14-16]. The device we present consists of the 
fullcrcnc molecule Ceo fixed in the middle of two gold electrodes. Here, we consider two distances between the 
electrodes, 11. tA and 13.7 A. This correspondes to the minimum distance between an Au atom and C atom of 2.3 A 
and 3.3 A, respectively. Each electrode consists of repeated unit cells with nine Au atoms in the (100) direction and 
extended to infinity. Note that we do not expect that the Ceo molecule can be found in the Coulomb blockade regime 
in such junctions if the electrode separation is small, say 9 — 13 A. The number of electrons in the molecule is not a 
good quantum number if the interaction with electrodes is strong. 
The current- voltage (I-V) characteristics is calculated as follows: 

1=^1 H) [fi{E) - fr{E)] dE , (6) 

where T{E, Vb) is the transmission coefficient at the energy E and bias voltage 14- The Fermi distribution functions 
/; (E) limit the integration in equation (6) to the small energy range eVb at the Fermi level of the electrodes (in the 

zero temperatiuc limit). Transmission coefficient T{E, Vb) is calculated by use of the Green's functions [17] which we 
have already obtained as a result of the NEGF-DFT self-consistent iterations, so we write 

T{E, Vb) = TV [Ti{E) G«(iJ) TriE) G« t(£;) ] . (7) 

Here, the matrix F; (E) is related to the line- widths resulting from the coupling of the scattering region to the left 
(right) lead. It is evaluated in terms of the corresponding self-energy matrices (() as follows: 

r;(,)(£) = i(sf(,)(£;)-i]f(t)(£)). (8) 

The retarded Green's function of the scattering region G^ {E) is given as the following matrix: 

G''{E) = {ES-F[p]-^\E))-\ (9) 

where F[p] is the effective Hamiltonian of the scattering region and — -I- is the total self-energy matrix. 
The density of state function is calculated from the density of state matrix A{E) = i {G^{E) — G^^E)) as follows: 

DOS{E) = Tr {A{E) S) /2'k. 

For the Au-Ceo-Au junctions at Vb = 0, the density of states near the Fermi energy is presented on figures 2 and 
3 (the Fermi level is set to Ep = in this article). Note that at positive gate voltages, Vg > 0, the energy levels 
originated from the LUMO of Ceo are partially occupied so they lie near E = and Q > (in the units of e). 
At negative gate voltages, Vg < V* < 0, the energy levels originated from the heighest occupied molecular orbital 
(HOMO) of Ceo begin to be depopulated, so the energy levels originated from the LUMO of Ceo are shifted far to 
the right from E = and Q < (in the units of e). 
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FIG. 2. Density of states of the Au-Ceo-Au junction and excess charge of the Ceo molecule are presented for the electrode 
separation of 11.7 A and VJ, — 0. At the zero gate, the excess charge is « 0.7 e. Note that the Ceo molecule is not in the 
Coulomb blockade regime (the broadening is large). 
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FIG. 3. Density of states of the Au-Ceo-Au junction and excess charge of the Ceo molecule are presented for the electrode 
separation of 13.7 A and VJ, = 0. At the zero gate, the excess charge is « 0.5 e. Note that the Ceo molecule is not in the 
Coulomb blockade regime yet. 
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FIG. 4. Transmission coefficient T{E, 14) versus electron energy E is presented for the electrode separations of 11.7 A and 
three bias voltages Vj,, Vg — 0. The sharp transmission peak, which results from a resonances of the fuUerene + contacts 
Hamiltonian, dominates the conductance as well as the current-triggered dynamics. Positions of the Renormalized Molecular 
Levels (RML) [15] of Ceo are depicted over the peaks of T{E, Vb), and these levels are the LUMO-derived states. The inset 
shows the calculated I-V curve from which the metallic behaviour of the junction is evident. 

The transmission coefficient obtained from equation (7) is presented on figure 4 as a function of tlie electron energy 
E for three different bias voltages and tlie zero gate potential. The sharp peak in T{E, Vb) is the result of a resonance 
transmission through the LUMO of the Cgo- At Vt = 0, this resonance lies just above the Fermi level of the leads, 
see figures 2 and 3. As the bias voltage is applied, the peak position shifts toward higher energies, ff the system is 
absolutely symmetric and Vb ^ 0, we would expect the same voltage drop at the two electrode-molecule junctions. 
In our system, the atomic structures of Ceo facing the left and right electrodes are not the same, and this breaks 
the left-right symmetry. We therefore observe an asymmetric voltage drop at the two contacts. Consequently, the 
transmission peak position shifts by about ~ O.lSeV when Vb is increased by 0.2V, i.e., the apphed voltage drops 
more on one side than on the other side of the Ceo- We also calculated the I-V curve for this device by use of equation 
(6); it is shown in the inset of figure 4. 
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FIG. 5. Equilibrium conductance G is plotted as a function of the gate volage Vg for the electrode separations of 11.7 A. 
Note that G ^ Go at Vg = 0. 
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FIG. 6. Equilibrium conductance G is plotted as a function of the gate volage Vg for the electrode separations of 13.7 A. 
Note that G ^ 0.25 Go at Vg = 0. 

We can easily calculate the equilibrium conductance G of the junction, G = dl/dV. In the linear regime, —0.25 V < 
Vb < 0.25 V, see figure 4, we can use the following formula: 

G = GoT{E^Ep,Vb = 0) (10) 

(valid in the zero temperature limit). Here, Gq = 2e^//i is the quantum of conductance. As the positive gate potential 
changes the occupancy of the LUMO-derived energy levels, we can expect the strong dependence of G{Vg) for different 
separations of the leads. These results are presented in figures 5 and 6. 

Similar metallic behaviour has been reported before [14] for Al-Ceo-Al junctions at the electrode separation of 9.3 A 
(it correspondes to the minimum distance between an Al atom and C atom of » 1.1 A). Note that the work function 
of Al is « 4.19 eV whereas the work function of Au is « 5.3 eV. At Vt = Vg = 0, we have Q « 3e, G « 2.2 60 for 
the Al-Ceo-Al junction [14] and Q « 0.7e, G « Go for the Au-Ceo-Au one (the electrode separation of 11.7 A). In 
addition, for the Au-Ceo-Au device, the current in the linear (metallic) regime is (roughly) one order of value smaler 
than the current calculated for the Al-Cgo-Al one. 

Here, the interesting physics is that despite the large HOMO-LUMO gap of a free Geo (— 1.8 eV [19]), we predict 
metallic transport characteristics for the Au-Ceo-Au devices, and they are qualitatively the same as the properties of 
the Al-Ceo-Al devices. We found that if the Geo molecule is relatively well bonded to metallic leads, there is a strong 
charge transfer from the leads to the Ceo cage and it partially fills the (empty) LUMO of Ceo- In other words, charge 
transfer from the electrodes to Ceo aligns the LUMO to the Fermi energy of the leads. As a result, we obtain a large 
resonance conductance and the metallic I-V curve. 



IV. SUMMARY 

We have shown that the NEGF-DFT formalism is a powerful technique for modeling charge transport properties of 
molecular electronic devices. The novelty of this technique is in constructing electron charge density via nonequilibrium 
Green's functions and the very effective screening approximation. Since the entire algorithm is based on evaluating 
Green's function, the technique is intrinsically 0{N) [20]: this is because atoms far away from each other do not 
overlap so that the Hamiltonian matrix is block-diagonal. This has tremendous computational advantage which we 
demonstrated by calculating the nonlinear I-V curve of the Au-Cgo-Au molecular tunnel junction. 
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